function varargout = interpdx(x_data, y_data, x_interp, num_deriv);

% INTERPDX - interpolate F(X) and calculates dF(X)/dX
%
%    YI = INTERPFX(X, Y, XI) estimates the values YI at the corresponding
%    values of XI from the vectors X and Y using spline interpolation.
%
%    [YI, DYI] = INTERPFX(X, Y, XI, 1) also returns the first derivative
%    of DYI/DXI.
%
%    [YI, DYI, D2YI ... DNYI] = INTERPFX(X, Y, XI, N) calculates all
%    derivatives of DYI/DXI up to the Nth derivative.
%
% Note that the spline interpolation method is selected because it is the
% only method that results in continuity of derivatives.

% By:   S.C. Molitor (smolitor@med.unc.edu)
% Date: August 4, 2000

% initialize outputs

for i = 1 : nargout
   varargout{i} = [];
end

% validate input arguments

if ((nargin < 3) | (nargin > 4))
   msgbox('Invalid number of arguments', 'INTERPDX Error', 'warn');
   return
elseif (~isnumeric(x_data) | isempty(x_data))
   msgbox('X must be a numeric array', 'INTERPDX Error', 'warn');
   return
elseif (~isnumeric(y_data) | isempty(y_data))
   msgbox('Y must be a numeric array', 'INTERPDX Error', 'warn');
   return
elseif (length(x_data) ~= length(y_data))
   msgbox('X & Y must have the same length', 'INTERPDX Error', 'warn');
   return
elseif (~isnumeric(x_interp) | isempty(x_interp))
   msgbox('XI must be a numeric array', 'INTERPDX Error', 'warn');
   return
elseif (nargin == 3)
   num_deriv = 0;
elseif (~isscalar(num_deriv))
   msgbox('N must be a scalar', 'INTERPDX Error', 'warn');
   return
elseif (nargout ~= num_deriv + 1)
   msgbox('Number of outputs must equal N + 1', 'INTERPDX Error', 'warn');
   return
end

% make sure X is monotonic
% determine if X is equally spaced for fast methods

diff_x = diff(x_data);
if (all(diff_x < 0) | all(diff_x > 0))
   if (min(diff_x) == max(diff_x))
      interp_method = '*spline';
   else
      interp_method = 'spline';
   end
else
   msgbox('X must be monotonic for interpolation', 'INTERPDX Error', 'warn');
   return
end

% interpolate values of YI

y_interp = interp1(x_data, y_data, x_interp, interp_method);
varargout{1} = y_interp;

% make sure XI is monotonic
% determine if XI is equally spaced for fast methods

diff_x = diff(x_interp);
if (all(diff_x < 0) | all(diff_x > 0))
   if (min(diff_x) == max(diff_x))
      interp_method = '*spline';
   else
      interp_method = 'spline';
   end
else
   msgbox('XI must be monotonic for derivative calculations', 'INTERPDX Error', 'warn');
   return
end

% iterative loop to calculate derivatives
% discrete derivative is secant between two pts
% shift x values to reflect this & reinterpolate
% endpoint is out of range & assigned NaN

x_deriv = (x_interp(1 : end - 1) + x_interp(2 : end))/2;
for i = 1 : num_deriv
   dydx = diff(y_interp) ./ diff_x;
   y_interp = interp1(x_deriv, dydx, x_interp, interp_method);
   varargout{i + 1} = y_interp;
end
return
